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Abstract 

We present a new particle-merging algorithm for the particle-in-cell method. 
Based on the concept of the Voronoi diagram, the algorithm partitions the 
phase space into smaller subsets, which consist of only particles that are in 
close proximity in the phase space to each other. We show the performance of 
our algorithm in the case of the two-stream instability and the magnetic shower. 
Keywords: particle merging, PIC, Voronoi, clustering, two-stream instability, 
magnetic shower, QED cascade 


1. Introduction 

For more than 60 years the particle-in-cell (PIC) technique [T] has been 
used to simulate a wide variety of physical problems, ranging from electrical 
discharge to particle acceleration. However, in several scenarios - in particular 
5 field ionisation or QED cascades - the number of particles in the simulation box 
grows exponentially. Due to an overwhelming number of particles, the associ¬ 
ated memory required can easily exceed that available on even high performance 
computers and as a consequence the computational performance drops drasti¬ 
cally. 

10 In these situations, a particle merging algorithm (PMA) has to be imple¬ 
mented. The main goal of a PMA is to reduce the number of particles in a 
simulation box while keeping the physical properties of the system intact after 
a merging event. A straightforward PMA is to randomly pick a pair of particles 
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and then merge, see for example [2]- Since it merges with no guidance, the 
15 method is not able to preserve the phase space distribution, and so the physical 
picture is likely to be distorted after merging. The problem is that it fails to 
incorporate the notion of proximity in the phase space, i.e. how similar particles 
are, into its framework. In the scope of this paper, we call this PMA the blind 
method. 

20 Lapenta already proposed a scheme for merging particles (called “particle 
coalescence”) in |3] and |1]. In this method, particles are first sorted into two 
bins. Then the binning process continues until the number of particles per bin is 
small enough for the pairwise comparison. This type of PMA was then refined 
and improved by Teunissen and Ebert in which the k-d tree method was 
25 employed to search for the nearest neighbour. Recently, a similar approach was 
also proposed by Vranic et al. [6], where the momentum space is divided into 
smaller subcells for sorting particles. 

We design our PMA from a different point of view, in which the algorithm 
not only merges particles which are close in the phase space but also offers 
30 users a direct control over errors introduced by a merging event. The notion 
of proximity in our algorithm is developed through the concept of the Voronoi 
diagram [3, thus the name Voronoi PMA. As shown later, the quantification 
of the error is realised through the coefficients of variation. The algorithm is 
successfully implemented into the framework of the VLPL (Virtual Laser Plasma 
35 Laboratory) code [5]. 

The paper is organised as follows: in sectionj^ we briefly introduce the defini¬ 
tion and some examples of the Voronoi diagram; the comprehensive description 
of our PMA is revealed in sectionin section]^ we test the performance of our 
merging algorithm with three cases: the counter-propagating plasma blocks, the 
40 two-stream instability, and the magnetic shower simulations; finally, we sum¬ 
marise the paper in section 
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2. Voronoi diagram 


For any given set of n sites, S = {si,S 2 , in the real d-space the 

Voronoi cell Vk associated with the site Sk is a set of points in such that the 
distance from those points to Sk is not greater than the distance to any other 
site Sj {j ^ fc) in S' [5]. 

Vfc = {x e K''* I Vj : dist(a;, Sk) < dist(a;, Sj)} for 1 < i,j < n. (1) 

Here, dist(a;, y) denotes the metric function of the distance in The Voronoi 
diagram was first developed, though informally, in 1644 by Descartes. In 1908, 
45 the Russian-Ukrainian mathematician G. F. Voronoi formally defined and stud¬ 
ied the general case [7j. The concept is used in many contemporary research 
fields, such as geophysics, meteorology, and condensed matter physics. 

Original Space Euclidean Measure Chebyshev Measure 



Figure 1: The Voronoi diagrams with different metric functions. Each Voronoi 
region is painted with a distinct colour. The black star in each region is the 
Voronoi centroid. 

Observing eq. 0 , we see that the metric function dist(a;, y) plays a vital role 
in the formation of the Voronoi diagram. Different metrics will result in different 
Voronoi diagrams. Moreover, in our case, different metrics also require different 
implementations of the algorithm (see section]^ for more detail). Fig. shows 
the Voronoi diagram of a random distribution with Euclidean and Chebyshev 
measures. Given two vectors p and q, the Euclidean distance is 

dist(p, q) = (2) 
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while the Chebyshev distance is given by 


dist(p,q) = max|gi - Pi\. 


( 3 ) 


3. Algorithm 

Input: The algorithm requires two user inputs, Tx and Tp, which are the 
tolerances for position and momentum. These parameters are employed as the 
stopping condition and appear at step 3. A merging event will take place in a 
simulation cell if the particle number N of that cell is greater than the minimum 
particle number Amin- 

Step 1: For every simulation cell, collect all particles (weight Wi, position x^, 
and momentum pi) in that cell into a set Vo- This set Vo is the hrst Voronoi cell, 
which covers the entire phase space of a simulation cell. We then calculate the 
statistical average in the phase space of this set of particles Vo by the following 
formulae: 


Wo = ^ 
ieVo 

^0 — -, 

E^evo 

p _ E»gVo 

-^0 


( 4 ) 

( 5 ) 

( 6 ) 


The point (Xq, Po) with weight Wq is the centroid of the first Voronoi cell Vo- 
From now on, quantities of a Voronoi centroid are denoted by the capital letters. 

Step 2: We calculate the standard deviation of each dimension I in the 
phase space with respect to the current Voronoi centroid: 



( 7 ) 

( 8 ) 


We compute the coefficient of variation (CV) A for each dimension. The CVs 



for spatial and momentum dimensions are defined as 
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( 9 ) 

( 10 ) 


For the spatial dimensions, due to the symmetry in space the CV Axo is defined 
as the ratio between the standard deviation and the length Lxo of the first 
Voronoi cell Vo- On the other hand, since there is no such symmetry in the 
momentum space, the CV Ap^ is obtained from dividing the standard deviation 
60 by the mean value. As the CVs are dimensionless numbers we can treat the data 
obtained from the position and momentum spaces equally (see step 4 below). 
In our algorithm, the CVs represent the accuracy of the merging scheme, with 
smaller CVs resulting in smaller errors due to merging. 

Step 3: We compare the recently obtained CVs Axq and Ap^ with their 
65 corresponding tolerances Tx and Tp. If a Voronoi cell has all six CVs less than 
or equal to the tolerances, the algorithm will mark that cell finished and stop 
dividing it. On the other hand, as long as there is at least one component 
whose CV does not satisfy the aforementioned requirement, the algorithm will 
keep going to the next step. 

Step 4: We consider the individual components of Axq and Ap(,, that is 
{Aj;, Aj,, Az, Ap ^, Ap ^, Ap^}, and find the axis k which has the largest deviation. 


fc = maxA;, with I e {x,y,z,px,py,pz}- (11) 

70 Step 5 : Make a hyperplane cut through the the Voronoi centroid perpen¬ 
dicular to the axis k. Denote q and Q the dynamic variables of the particles and 
of the centre, respectively, on the axis k. The hyperplane cut divides the set Vo 
into two new independent subsets Vi and V 2 , whose new centroids are given by 
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Vi = {i G Vo : qi < Q} 

V 2 = {i GVo ■■ Qi 

= EzgVi 

= E7GV2 

X, 






E.gvi-. 



75 Step 6: Sort the particles into their corresponding new sets. Repeat steps 
2-6 for the new sets Vi and V 2 until the stopping condition is satisfied. 

Step 7: If the stopping condition is met for all Voronoi cells, the algorithm 
removes all particles from the simulation cell and replaces them with the Voronoi 
centroids as the merged particles. The algorithm ends here. 

80 We have several remarks on our algorithm: 

• Our Voronoi PMA is inspired by Schreiber’s adaptive k-means clustering 
algorithm used in Computational Geometry [10]. 

• In step 1, we state that the merging process is carried out cell by cell. 
However, the algorithm can be adjusted such that the first Voronoi cell 

85 Vo contains all particles from the simulation box and starts merging from 

there. The rest of the algorithm is kept intact. However, it is likely that 
the global merging approach violates the local charge conservation. In 
this case, one must take into account a correction scheme in order to com¬ 
pensate for the error caused by merging events. Which implementation is 

90 used depends strongly on the user preference or the code framework. We 

adhere to the cell-by-cell implementation as it is readily parallelised. 

• The distance measure used here (see eq. 0. step 4) can be considered as 
a Chebyshev-like distance, since eq. is not guaranteed for every particle 
and phase space dimension. We have chosen this measure instead of a 

95 more obvious candidate, the Euclidean measure, for the following reasons: 

1. The simplest implementation of the Euclidean measure requires the 
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seeding of Voronoi centroids at the beginning of the algorithm. More¬ 
over, the number of Voronoi centroids is kept constant throughout 
the algorithm. This limitation not only reduces greatly the flexibility 
100 of the algorithm but also cannot fit well to the dynamic situation of 

a physical problem uni. Conversely, the Chebyshev measure requires 
no seeding and suits perfectly the divide-and-sort scheme, which is 
applied here. 

2. In the author states a rule of thumb that for a given dataset of 

105 N points, the number of centroids is set to fc « •y/iV/2. Again, the 

number of Voronoi centroids cannot be changed once the algorithm 
starts. As such, we do not follow this rule. 

3. In order to use the Euclidean measure without a fixed number of 
centroids, we would have to solve the problem of an unknown number 

110 of clusters in a dataset. This can be done through the Bayesian 

information criterion m or the removing centroids method |13j . The 
former approach is difficult to implement, while the latter tends to 
be computationally intensive. 

• In the momentum space, the Voronoi PMA groups particles by taking into 
115 account both the direction and the magnitude of particles’ momenta. Due 

to the difference in the direction, it might occur that the energy is lost 
after a merging event. The relative error in the total energy is observed in 
Fig. [6] for the two-stream instability and Fig. [^for the magnetic shower 
below. These graphs show that the loss in energy per merging event is 
120 extremely small. However, the merging quality can be further improved by 

introducing a mechanism to conserve energy perfectly and directly. One 
can consider the Langdon-Marder corrector-scheme mmm or follow 
the proposal to merge into two particles [5]. We also make a side remark 
that the Langdon-Marder scheme becomes obligatory in case users want 
125 to implement the algorithm through the global merging approach. 
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4. Simulation 


Having presented the algorithm, we proceed to test its performance. To this 
end, we consider three situations: counter-propagating plasma blocks, the two- 
stream instability HZl [IE], and the magnetic shower produced by an energetic 
130 particle entering a strong magnetic field Hu- 

Before going further, we briefly describe the implementation of the blind 
method used here for comparison. We define the parameter a as the merging 
fraction. A merging event will take place in a simulation cell if the number 
of particles N of that cell satisfies the condition N > ceil(aA^). Then, the 
135 blind method merges particles in the current cell until the number of particles 
after merging is at maximum ceil(aA^). This implementation allows the blind 
method produces the same number of particles as in the Voronoi PMA for fair 
comparison. 

4-1- Counter-propagating Plasma Blocks 

140 The counter-propagating plasma blocks simulation is a simple test, in which 
two blocks of non-interacting particles with uniform density distribution prop¬ 
agate and then overlap each other. These blocks have the same momentum 
magnitude but opposite propagation directions (see Fig. [^. With no merging, 
there is no change to the system apart from the translation in x-direction af- 
145 ter the blocks pass through each other. By using this test we can easily spot 
whether a given PMA preserves the phase space distributions since there is a 
duration when the blocks overlap. If a merging method does not preserve, two 
or more particles from the different distributions might be merged together. 
Here, we compare the performance of the Voronoi PMA and the blind method. 
150 The merging period Tmrg = 2At, with At is the time step, is applied for both 
methods. For the Voronoi PMA, the tolerances are Tx = 0.4 and Tp = 0.01. 
For the blind method, we deliberately choose the parameter a so as to give a 
similar final number of particles as in the Voronoi PMA. 

We look at the number of PIC particles appearing in the simulation (see 
Fig. §1). Starting with 12000 particles, the blind method merges into 4200 
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Figure 2: The initial configuration for the counter-propagating plasma blocks 
simulation: two blocks have the same momentum magnitude but opposite prop¬ 
agation directions. The merging event will commence when two blocks start 
overlapping each other since the particle number exceeds the threshold. For 
the Voronoi PMA, the threshold is Nmin = 15, and ceil(aiV) for the blind 
method. A good merging algorithm will leave behind no change in the phase 
space distribution apart from the translation in the x-direction. 


particles at the end of the simulation, while the Voronoi PMA finishes the 
task with approximately 3800 particles. The numbers of particles produced 
by two methods are approximately equivalent. Fig. shows the phase space 

distributions at the end of the simulation and figs. (b,c, and d) show the 

160 histogram. For the blind method, we see that after the blocks have passed 
through each other, there are many particles left behind between the two blocks. 
The momentum space plot and the histogram shows that these particles have 
zero momentum. The blind method also produces many particles with momenta 
not equal to the original magnitude (150toc). As a consequence, the particle 
165 distributions are smeared and the conservation of energy is violated. Conversely, 
the Voronoi PMA accurately preserves the phase space distributions, returning 
the same result as for the case with no merging. For this test, we see that 
despite the fact that it finishes the simulation with fewer particles than the 
blind method, the Voronoi PMA accurately preserves the particle distributions. 
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Figure 3: The number of PIC particles during the simulation (fig. a) and 
the histograms for no-merge, the Voronoi PMA, and the blind method (figs, 
b, c, and d respectively) in the counter-propagating plasma blocks simulation. 
Despite merging into a similar number of particles, the Voronoi PMA does not 
distort the momentum distribution. 

170 while the blind method does not. 

4 . 2 . Two-stream instability 

The two-stream instability consists of two identical particle beams stream¬ 
ing through each other. These beams propagate in the opposite directions and 
a small perturbation in the charge density can change the electric field, which 
175 in turn causes further perturbation in the density distributions. This type of 
simulation makes an illustrative example of how the algorithm manage merging 
particles in a dynamic evolution of the phase space. The configuration for the 
two-stream instability is listed in table At the beginning of the simulation, we 
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Figure 4: The phase space distributions (first row x/y, second row xjpx) at 
the end of the counter-propagating blocks simulation. The blind method leaves 
behind many particles that have zero momentum. Meanwhile, the Voronoi PMA 
reproduces the result obtained with no merging. 


create two electron beams with the same initial Lorentz factor 7 = 1 but oppo- 
180 site propagation directions. Each beam has 16 x 10"^ particles and is neutralised 
by the background charge density. Purposefully, the merging algorithms are 
only enabled after time t = 5Ao/c, when the instability can be visibly observed. 
The merging fraction for the blind method is chosen to be a = 0.965, such that 
we can have a fair comparison between two algorithms. 


Wavelength 
Simulation box 
Grid steps 
Time step 

Electron initial Lorentz factor 
Number of CPUs 
Merging period 


Aq = 800 nm 

3 . 2 A 0 X I.OAq 
O.OIAq X O.IAq 
A t = 0.005Ao/c 

7 = 1.0 

8 X 1 
50At 
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Merging start 

5Ao/c 

The minimum particle number per cell 

200 

(for Voronoi PMA) 


Tolerances (for Voronoi PMA) 

Tx = 0.8 and Tp = 0.15 

Merging fraction (for the blind method) 

a = 0.965 


Table 1: The configuration for the two-stream instability simulation. 


185 

The phase space distribution {xjpx) for the two-stream instability is shown in 
Fig. 0 at different time stamps. Similarly to the counter-propagating plasma 
blocks, the blind method (the last column) produces many particles with mo¬ 
menta approximately equal to zero, which do not appear in the original simu- 
190 lation (the first column). This early distortion in the phase space distribution 
leads to a different instability growth at later time. On the other hand, the 
Voronoi PMA (the second column) retains the phase space distribution through¬ 
out the simulation. Moreover, in contrast to the smooth pictures obtained with¬ 
out merging, the outcomes of the two algorithms appear grainier, since there are 
195 lesser particles in the phase space due to merging events. Fig. [^shows the num¬ 
ber of electrons and the relative error in the total energy 5E = abs(iilo — E)/Eq. 
Here, Eq is the energy of the system without merging. Observing Fig. [^, we 
see that when the merging event is enabled (at t = 5Ao/c), there is a steep fall in 
the number of particles for the Voronoi PMA (the red line), falling from 32 x 10"^ 
200 to 20 X 10^ particles. This abrupt drop is then followed by a short decline to 
17.8 X 10^ particles. At around t = SXq/c, there is almost no merging event till 
the end of the simulation, since the number of particles per cell is already below 
the threshold. On the contrary, the blind method (the green line) exhibits a 
steady decline in the number of particles , reducing to 16.8 x 10^ particles at 
205 the end of the simulation. From Fig. |6]3, wee see that the total energy relative 
error is rising up to 0.1 for the blind method, while the Voronoi PMA reaches 
a peak at 5E = 0.006 during the simulation. 
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Figure 5: The phase space distributions (x/px) for the two-stream instability 
simulation at different time stamps. The first column shows the original simula¬ 
tion with 32 X 10^ particles. The second and third columns show the simulation 
with the blind and Voronoi merging method, respectively. While the outcome 
produced by the blind method looks different, the Voronoi PMA follows the 
evolution course as in the no-merging case. 
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Figure 6 : The number of PIC particles during the simulation (fig. a) and 
the relative error in the total energy due to merging events for the two-stream 
instability simulation. The Voronoi PMA reduces the number of particles from 
32 X 10^ to 17.8 X 10^ particles and stops merging from there, since the number 
of particles per cell is already below the threshold. The highest relative error in 
the total energy for the Voronoi PMA is 0.006, and the blind method 0.1. 

4-3. Magnetic Shower 
4-3.1. Introduction 

Consider an energetic particle propagating through a strong magnetic field. 
Due to the interaction with the field, the particle will emit hard photons on 
its course. In turn, these photons interact with the field and will decay into 
energetic electron-positron pairs. The cascade of particles develops quickly and 
an exponential growth of the number of particles is usually observed. This 
phenomena is called the magnetic shower. The occurrence of the magnetic 
shower requires both an intense field and high particle energies [19] [20]. This 
condition is quantified in the quantum parameter y [19j . which is defined as 

B , , 

X = l-^- ( 12 ) 

Bs 

Here, 7 is the particle’s Lorentz factor, B the magnetic field strength, and 
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the Schwinger field Bs = 4.41 x 10^^ G. The pair production has sufficient 
probability to start the cascade process only when y > 0.1 [19]. The proba¬ 
bility rates for photon emission and pair production are expressed in intricate 
expressions (see eq. ( 2 ) and (3) in ref. |1T| and the description therein). The 
computation usually requires solving the double integral of the Airy function. 
Thus, the task involves a significant computational overhead. However, under 
the assumption that the dimensionless field amplitude Oq ^ 1 , the field can be 
regarded as being constant during the decay processes. Additionally, if both 
conditions y ^ B/Bs and B ^ Bs are satisfied, we can utilise the theory 
of quantum processes under a constant cross field given in |52| [13|. Accord¬ 
ing to this theory, the probability rates for the photon emission Wem and pair 
production Wpair are 




a mc^ f^hx^ + lx + 1 f2x\ 
Sv^TT h'y Jq (l-fa:)3 V^yj 


(13) 


and 
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3(1 — x'^)k 


da;. 


(14) 


210 Here, a is the fine structure constant; ^ 2 / 3 (a;) is the modified Bessel function 
of the second kind; e is photon’s energy and k its quantum parameter. Our 
numerical model for the cascade process is based on the Monte Carlo method 

m m- 


The magnetic shower is an appropriate example since the number of particles 
215 can grow exponentially during the simulation and the particles’ energies can 
range from several to hundred MeVs. Thus, it is a good indicator of how a PMA 
copes with the dynamic development during the simulation while preserving the 
physical features of the system. 


4.3.2. Simulation 

220 The simulation configuration for the magnetic shower is listed in table 
We begin the simulation with 5 numerical electrons. For an electron with a 
Lorentz factor 7 = 5 x 10^ and a magnetic field eB/nieCU! = 500, the quantum 
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parameter is x « 150 ^ 1. Here, e is the elementary charge, rrie the electron 
mass, c the velocity of light, and w = 2'kcIXq. As before, we consider three 
225 cases: without merging, with the blind merging method, and with the Voronoi 
algorithm. As before we deliberately choose the merging fraction a such that 
the blind method and the Voronoi PMA result in the similar number of particles 
at the end of the simulation. 


Wavelength 

Ao = 800 nm 

Simulation box 

' 3.2Ao X 3.2Ao x 3.2Ao 

Grid steps 

0.04Ao X 0.04Ao x 0.04Ao 

Time step 

At = 0.005Ao/c 

Magnetic field strength 

H = 6.6 X 10i° G 

Electron initial Lorentz factor 

7 = 5 X 10^ 

Number of CPUs 

5x5x5 

Merging period 

2At 

The minimum particle number per cell 

10 

(for Voronoi PMA) 


Tolerances (for Voronoi PMA) 

Tx = 1.0 and Tp = 0.02 

Merging fraction (for the blind method) 

a = 0.88 


Table 2: The configuration for the magnetic shower simulation. 


230 The growth in particle number is shown in Fig. Without merging (blue), 
both electron and positron display exponential growth during the simulation. At 
the end of the simulation, a total number of 2.8 x 10® particles has been reached 
for each specie. Meanwhile, the photon specie grows from 0 to 7 x 10^ particles 
at the last frame. The blind method (green) results in 1.45 x 10® electrons and 
235 posittrons, 4000 photons. The Voronoi PMA (red) produces in total 1.35 x 10® 
electrons and positrons, and 8000 photons. That is, the number of particles in 
the box is reduced approximately 40 times by both methods. In order to verify 
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Figure 7: The number of particles in the magnetic shower simulation as a func¬ 
tion of time for electron, positron, and photon (from left to right). The first row 
shows the result from the simulation without merging, the second row shows 
the outcome by using the Voronoi PM A (red) and the blind method (green). 
The Voronoi PMA reduces the number of electrons (positrons) from 2.8 x 10® 
to 6.5 X 10^ particles. 

the validity of the simulation, we look at the total energy and the spectra of 
the particles. Figs. and illustrate the evolution of the particle energies and 
240 their spectra at the end of the simulation. For the blind method (solid, green 
line in Fig. [^, we see a gradual decrease in the total energy of electrons and 
positrons around the point when the photon energy is reaching its peak. This 
strongly affects the spectrum of every specie in the simulation box (see Figs, 
g, h, i): the distinct peak electrons and positrons is not observed. On the 
245 other hand, with a careful approach the Voronoi PMA (short dash, black line) 
overlaps the case with no merging (long dash, light blue )in Fig. showing 
that it preserves the physical behaviour in the total energy, with the decrease 
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in electron energy, increase in positron energy, and the sharp rise followed by a 
decrease in photon energy. Moreover, the Voronoi PMA accurately reproduces 
250 the spectra obtained with no merging (see Figs, [^d, e, and f). Originally, the 
simulation with no merging takes approximately 2 hours (7265 seconds). With 
the same settings, the Voronoi PMA completes roughly in 20 minutes (1172 
seconds) and the blind method takes about 24 minutes (1440 seconds). 



Figure 8: The total energy evolution of electron, positron, and photon during 
the magnetic shower simulation for three merging cases: no-merge case (long 
dash, light blue); Voronoi (short dash, black), and blind (solid, green). Unlike 
the blind method, the Voronoi PMA reproduces the results obtained from the 
original simulation. 


Finally, we perform a parameter scan on the tolerances Tx and Tp in order 
255 to observe the growth of particles and the accumulation of error due to merging. 


Fig. 10 shows the number of electrons and the relative error 5E = {Eq — E)/Eq 
during the simulation and Fig. displays the total computation time with 
various tolerance settings. Here, Eq is the energy of the system without merging. 
The most accurate simulation is achieved with Tx = 0.5 and Tp = 0.005. With 
260 this setting, the simulation takes roughly 40 minutes to complete and the total 
energy loss is around 0.05 MeV {6E « 1 x 10“^). We observe that the growth 
is also exponential and the number of electrons has reached 2.4 x 10® particles 
at the end of the simulation. When we loosen the tolerances, more particles 
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Figure 9: The spectra for the electron, positron, and photon species in the 
magnetic shower simulation at time t = 98At for three merging cases: no merge 
(blue), the Voronoi PMA (green), and the blind method(red). The spectra of 
particles are accurately reproduced by using the Voronoi PMA. Meanwhile, with 
the blind method, the distinct peak for electrons and positrons is not observed. 

are merged together. As a result, the growth rate becomes more linear but the 
265 energy loss develops speedily. In our test, the extreme case with Tx = 1.0 and 
Tp = 0.03 produces 7.7 x 10^ electrons and positrons, 5 x 10^ photons, and 
takes 14 minutes to hnish. However, in this case, it accumulates 20 MeV total 
energy loss {6E « 3.9 x 10“®). Although the loss is extremely small, we notice 
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the double in magnitude just by increasing from Tp = 0.025 to Tp = 0.03. We 
270 also observe that, the purple line (Tx = 0.5 and Tp = 0.02) completely overlaps 
the dark blue line (Tx = 1.0 and Tp = 0.02), showing that the tolerance Tp 
is more sensitive than Tx. Since, in a given cell, the particle momenta may 
vary significantly, an accurate simulation requires small Tp. We recommend 
Tp = 0.05 and Tx = 1.0 as a threshold for this type of simulation. 



Figure 10: The number of electron and the relative error in the total energy due 
to merging with various tolerances [Tx,Tp] settings for the magnetic shower 
simulation. With relaxed tolerances, the growth of particle number becomes 
linear but the error also accumulates faster. When stricter tolerances are used, 
the growth resumes the exponential behaviour while the error develops with a 
slower rate. We also observe that, the purple line (Tx = 0.5 and Tp = 0.02) 
completely overlaps the dark blue line (Tx = 1.0 and Tp = 0.02), showing that 
the algorithm is always more sensitive towards the momentum space. 
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Figure 11: The total computation time for the magnetic shower simulation. 
The last six columns show the simulations with the Voronoi PMA for various 
tolerances [TxjTp] settings. 

275 5. Summary 

In this paper, we present the Voronoi particle merging algorithm for PIC 
codes. The phase space of a simulation cell is partitioned, as in the Voronoi 
diagram, into smaller subsets, which only consist of particles that are close to 
each other. The quality of a merging event is ensured by two user inputs, the 
280 tolerances on position 7x and momentum Tp. The tolerances act as the balance 
between the speed-up and the accuracy of the simulation. Stricter tolerances 
mean smaller error but without much in the speed-up. On the other hand, re¬ 
laxed tolerances result in more merged particles and thus the computation time 
decreases but the error will accumulate faster. Making a right combination for 
285 the tolerance pair for a certain simulation requires prior knowledge of particles’ 
behaviour. If a simulation involves particles which spread out in a large range 
in the momentum space, we suggest keeping the Tp lower than 0.02. Other¬ 
wise, this value can be relaxed. On the other hand, since it relates to particles’ 
relative position in a cell, Tx can be chosen up to 1.0. 
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290 We have tested the performance of our algorithm with three tests: the 
counter-propagating plasma blocks, the two-stream instability and magnetic 
shower simulations. In all cases, we observe that the conservation of momen¬ 
tum is perfectly held and the conservation of energy is maintained extremely 
well, with only small margin of error. The two-stream instability shows that the 

295 Voronoi PMA preserves the phase space evolution and the total energy error 
in this case is of the order of 10“^. In the magnetic shower simulation, the 
total energy error is of the order of 10“® with a speed-up by a factor of 6 and 
the spectra of particles are also comparable very well to those obtained with no 
merging. 
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